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Abstract. 

A 3D model intermediate between cellular automata (CA) models and the reduced magnetohydrodynamic 
(RMHD) equations is presented to simulate solar impulsive events generated along a coronal magnetic loop. 
The model consists of a set of planes distributed along a magnetic loop between which the information propa- 
gates through Alfven waves. Statistical properties in terms of power-laws for energies and durations of dissipative 
events are obtained, and their agreement with X-ray and UV flares observations is discussed. The existence of 
observational biases is also discussed. 
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1. Introduction 

One of the main unresolved problems in solar physics con- 
cerns the mechanism by which the solar atmosphere is 
heated from several thousands degrees in the photosphere 
to millions in the corona. It is now commonly accepted 
that the ultimate source of energy lies in the convective 
motions in and below the photosphere, and that a reliable 
model of coronal heating has to deal with the transfer, 
the storage and finally the release of this energy into the 
solar corona. Several conceptual models have been pro- 
posed, such as Alfven waves, electric curr ents and MHD 
turbulence - see for a review IZirkerl l)l993|) - where in all 
cases the magnetic field plays a key role in the dynamics. 

On the other hand, it is also well established 
that impulsive events (e.g. solar flares, X-ray 
bright points) are distributed in the corona over a 
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large events seem to be made up of the sup erposition of a 
myriad of smaller unresolved events. It was IParkerl (^988) 
who suggested that the active X-ray and UV corona is 
composed of a swarm of localized impulsive bursts of 
energy called nanoflares (or picoflares). Much theoretical 
work has been done to investigate Parker's conjecture and 
more generally the statistical nature of the solar coronal 



heating. They mainly follow two complementary schools 
which refer on the one hand to the dynamics of complex 
systems and on the other hand to fluid mechanics. 

The statistical properties of flaring activity 
allow one to view the solar corona as a com- 
plex system which can be described with cellu- 



lar autom ata ( CA ) models llLu fc Hamilton! Il99~ ; 
huetall Il993l IVlahos etall ll 99.4 ICals^aTc . 
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CA models have become increasingly useful in the 
study of complex systems because they permit the 
study of an entire system without ignoring the effects 
of individual components of the system. There are 
many natural applications such as substorms in the 
magnetotail i Takaloetall Il999lh st ar fo rmation in 
spiral galaxies llLeieune fc Perdanel Il996ll or earth- 
quakes ( Carlson fc Langert Il989l) . A cellular automaton 
is based upon the idea of the locality of influence : a 
system is distributed in space, and nearby regions have 
more influence than those far a part (see for instance 
iMacKinnon fc Macphersonl l|l997|l for a study of a nonlo- 
cal communication). A grid of cells is used to represent 
the components of a system, and each cell is given a set 
of phcnomcnological rules concerning its surrounding 
neighbors. The system evolves over several iterations by 
allowing each cell to interact using the given rules. What 
makes CA so interesting and useful is that after many it- 
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erations they reveal complex structures and arrangements 
that form across great distances even though each cell 
only takes into accou nt local informatio n . Self-Organized 
Criticalitv (SOC1 JBak et all Il987t iKadanoff et all 
Il989t iHwa fc Kardarfll992l: ISornettel l2000j) refers to the 
spontaneous organization of such an externally driven 
system into a globally stationary state over many scales. 

On the other hand, we find fluid models which give the 
physical description that is missing in CA. Much work has 



ocen none on statistical soiar nares (rjonecoDC &z oucian 
199i IWalsh, Bell & Hood, 1995: Einaudietall [l996 


Galsgaard & Nordlundl Il996 


: iHendrix & Van Hovcn 


19961 iDmitruk et all Il998l 


Galticr & Pouquen. Il998l 



2000) but most of them suffer from the fact that sta- 



tistical simulations of flares studied in the context of 
forced resistive MHD equations are possible only at the 
cost of huge computational expenses. Nevertheless it 
has been possible to show important properties, e.g. 
that the dissipative events produced exhibit power-law 
distributions (for total energy, peak of luminosity and 
duration) in agreement with X-ray observations, but with 
generally a much smaller "inertial range" than the CA 
counterpart. 

A recent debate about the possible existence of 
sympathetic flaring, i.e. the correlation in t i me of 
two successive events jPearce. Rowe &; Yeund. Il993l 



Wheatland, Sturrock & McTiernan, 1998; Boffctta et al 



1999t IWheat andl l2000t iLenreti. Carbone fc Veltril l200lt 
Galtieii |2001 ). suggests the possibility to dismiss CA as 
a model of solar flares since standard CA models do not 
produce correlated events (non-Poissonian statistics such 
as power-law waiting time distributions). But in fact many 
CA mo dels exist in the liter a ture like nonconservative 
models l|Christensen fc Olamii Il992h which turn out to 
be able to generate the statistics expected for sympathetic 
flaring. But the question of the existence of sympathetic 
flaring in the corona has not yet found an answer since in 
particular there is still a debate about what we mean by 
event. 

The problem of coronal heating is intimately linked to 
the existence of nanoflares whose Probability Distribution 
Function (PDF) in energy is supposed to be a power-law 
steeper than that for regular flares. Let us assume that 
the PDF in energy E of events is distributed according 
to a power-law of index — f, i.e. Pr(E) oc E~**. It is then 
possible to show that ther e exists a critical slope of in- 
dex ( c = 2 i Hudsonl Il99lh . Indeed, the total energy re- 
leased in the corona by events between E m i n and E max 
is (E%& - E 2 -£)/{2 - which means that if C < 2 
the main contribution comes from high energy events, 
whereas if £ > 2 it comes from smaller events (the swarm 
of nanoflares). The average power dissipated in a large 
flare is of the same order of magnitude as the total aver- 
age power emitted by the corona, ~ 10 3 W • m -2 , which 
proves that regular flares can not account for coronal heat- 
ing since they are episodic events seen over and above the 
average coronal background. It is then natural to think 



that a swarm of very small and still unobservable events 
may dominate the heating process. One of the main chal- 
lenges of statistical flare models is to know whether or not 
it is possible to produce power-law distributions for any 
relevant quantity such as energy, luminosity or duration, 
and what the power-law indices are. 

The aim of this paper is to introduce a hybrid model 
for a solar magnetic loop which is somewhat intermedi- 
ate between CA models and full MHD or reduced MHD 
models. In this model, we will inject and store energy into 
a coronal loop (our numerical domain) via wave propa- 
gation from the photosphere (our numerical boundary). 
The trigger for an event is determined in a way analogous 
to conventional CA models, i.e. with a threshold in the 
current. However, during the subsequent event the cur- 
rent is dissipated and the magnetic field recomputed us- 
ing Maxwell's equations. Let us note that this is a minimal 
consistency requirement for the field evolution which is not 
always incorporated in CA models. In practice, the model 
allows current concentrations to form kinematically (ad- 
vection from the photosphere), but not dynamically (the 
nonlinear part of the Lorenz force, j x b). The model is 
non-trivial because of the threshold dynamics of the dissi- 
pation, which mimics the nonlinear terms, but the model 
is much simpler to integrate than the full MHD equations, 
therefore allowing a fast computation of statistics (events 
sizes and durations), and a comparison both with obser- 
vations and full numerical simulations. 

The paper is organized as follows. In Section 2 we give 
a detailed description of the CA model and show its basic 
behavior through some numerical experiments. In Section 
3 the results of a parametric study are given and discussed. 
In Section 4 we summarize the properties of the model, 
we present a comparison with observations, and we draw 
some conclusions. 

2. The model 

In the origin a l 3D lattice model developed by 
iLu fc Hamiltor] Jl99ll) the physical quantity defined 
on each lattice is the magnetic field. The system is 
driven from the outside by adding randomly in space 
a random magnetic field. The process continues until a 
reconnection instability criterion is satisfied at any point 
of the 3D lattice, i.e. until the magnetic gradient exceeds 
a critical value at this point. Then the magnetic field 
is redistributed (diffused) towards neighboring nodes 
with the possibility to transfer the instability as well. 
The redistribution process stops when the system is 
completely relaxed. Then another random amount of 
magnetic field is added to the system. An event called 
avalanche is associated to the rapid diffusion of the 
magnetic field. 

Subsequent models use the magnetic vector potential 
A rather than the magnetic field since the divergence- 
free condition for the magnetic field is then automati- 
cal ly satisfied. For example in the recent mode l developed 
by llsliker. Anastasiadis fc Vlahosl l)200fl l200lh where the 
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3D lattice represents an ensemble of magnetic loops, the 
knowledge of A allows to reconstruct the magnetic field 
topology and eventually the structure of the current den- 
sity. To do so they introduce the notion of derivative. The 
present simplified model belongs to this class of models 
but only one typical coronal loop will be considered and 
simulated. The detailed description of the model is now 
given. 



appears to be constant along the loop and is much smaller 
than their length. 



2.1. Basic idea: on-off mechanism and turbulence 

A possible reason for the behavior of the corona is that it 
lies in a turbulent state. A model of coronal loops should 
therefore allow for the effects of turbulent fluctuations. 
This is possible with CA models at a very superficial level 
through an on-off mechanism. The idea behind the thresh- 
old dynamics of our model, the on-off mechanism, is the 
following. The forcing due to the convective granules, al- 
though applied on a range of scales, has a typical length 
scale that is supposed to be far greater than the dissipative 
scale. The connection between the forcing and the dissipa- 
tive length scales is made through a turbulent mechanism. 
During the "off" phase, i. e. the loading phase, the plasma 
is in a laminar state where the dynamics is essentially gov- 
erned by the linear terms (and the loading). Because the 
system is driven slowly, parts of the system or even the en- 
tire system can stay in principle in this state for very long 
periods of time. When sufficient energy is accumulated in 
the loop some nonlinear instability appears which triggers 
the rapid generation of small scales. The inertial range 
of the turbulent energy spectrum extends to small scales 
and makes the link between the typical forcing length scale 
and the dissipative length scale. This "on" phase is there- 
fore characterized by a sudden increase of the dissipative 
terms in the RMHD equations (see section|^2J leading to 
a bursty event. Then the system returns immediately to 
an "off" phase. The nonlinearities of the RMHD equations 
are therefore included in the model in a very schematic 
way through a threshold dynamics only. 

2.2. Description of the model 

Geometry of the model. The model describes a coronal 
loop anchored in the photosphere whose footpoints are 
randomly moved. The presence of a strong axial magnetic 
field leads to essentially 2D dynamics, i.e. perpendicular 
to the mean magnetic fie ld, for which t he approximation 
of the RMHD equations (jStraussl Il976|) is a good model. 
As we can see on Fig.^ the 3D regular grid is made up of a 
set of planes distributed along the loop between which the 
information propagates through Alfven waves. Therefore 
each plane will evolve essentially independently from each 
other. Both boundary planes represent the photospheric 
footpoints, while the intermediate planes represent the 
loop itself, as if it were unbent. The curvature of the loop is 
not taken into account since the width of observed coronal 
loops (see e.g. observations from the TRACE instrument) 
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Fig. 1. Model of coronal loop used for the simulations: 
the loop is unbent into a box, whose both extreme planes 
represent the photosphere. Parallel planes represent slices 
orthogonal to the local direction of the loop. 



RMHD equations. Our aim is to compute the temporal 
evolution of the velocity field v and the magnetic field 
B (or b = Bj y/poHo if we consider only fields with the 
same physical dimension). We assume the presence of a 
strong and uniform axial magnetic field along the z-axis 
bo = boe z (= Bo/ y/pono) and we consider small pertur- 
bations to this field. We separate b—bo and v into parallel 
components (b z e z and v z e z ) and orthogonal components 
(b± and v±). With the following additional hypotheses, 
(i) the scales along the z-axis are larger than the scales in 
the orthogonal directions (gradients along the z axis are 
negligible), and (ii) the kinetic pressure is negligible com- 
pared to the magnetic pressure (i.e. the plasm a is cold , 
<C 1), we then obtain the RMHD equations ijStraussl . 

SI 

d t vj_ + (vj_ ■ V ± )t> ± = b a d z b ± + v A ± v ± (1) 

+(b ± ■ V±)b x - V(6i/2) , 

d t b x + • V_l)6_l = b d z v x + V (2) 

+(&±-V_l)«l, 

where v and r\ are respectively the kinematic viscosity and 
the magnetic resistivity. To each grid point are associated 
two scalar fields a s , with s = ±, from which the Elsasser 
fields are derived 

z s =v±+ sb± = V_l x a s e z . (3) 

All other fields (magnetic and velocity fields, current den- 
sity, vorticity. . . ) are derived from a s by analogy with the 
standard magnetic and kinematic equations. 

Initial state and boundaries. The fields a s are taken to be 
zero initially. Each cross-sectional plane along the loop is 
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periodic for the a s variables. All numerical computations 
for each plane are made in the Fourier space. We use in- 
verse Fast Fourier Transforms (FFT) when we temporarily 
need to know the values of a field at some positions in real 
space. 

Alfven wave propagation. In the right hand side of equa- 
tions (11121) . the first terms correspond to the Alfven waves 
propagation along the z-axis : a + propagates to the bot- 
tom of the simulation box (low values of z), while a~ 
propagates to the top. These propagations are modeled by 
specific cellular automaton rules used at each time step St 
of the simulation, corresponding to discretization of the 
Alfven waves terms of equations illlL'l) in the form : 

a s (z,t + 5t) = a s (z + sSz,t) (4) 

There is no loss of energy (or reflection) during the Alfven 
wave propagation since the density is assumed to be con- 
stant ; however we assume that there is a total reflection 
of the waves when they reach the photosphere, i.e. the two 
boundary planes of the simulation box. 

Loading. The action of the photospheric granules on the 
magnetic footpoints is modeled by a random increment 
S^ to the fields on both boundary planes. This incre- 
ment has random Fourier coefficients but has globally 
a power-law energy spectrum in ~ k~ a (the total in- 
tensity Pioad = J(fc<5 x f') 2 dfc is also a parameter of the 
simulation). Indeed, observational evidence suggests that 
the convective layer is in a turbulent state : the photo- 
spheric granules exhibit a turbulent power spectrum of 
velocity consistent with a Kolmogorov energy spectrum 
in fc -5//3 but only for a narrow i nertial range of wavenum- 
bers (I ~ 1/fc < 3 arcsecl ilRoudier fe Mulled. Il986t 
IChou et all Il99lt lEsoagnet et all Il993[) . We emphasize 
that there is no loading in the other parts of the loop : en- 
ergy is solely carried by Alfven waves. Furthermore, these 
waves reflect on both photospheric boundary planes. 

Dissipation criterion. We assume that dissipation occurs 
when an instability criterion is satisfied, which is the con- 
dition that the magnitude of the current density || J\\ ex- 
ceeds a critical value J c . As J z = Je z can be derived from 
the computed variables a s (J x and J y are negligible), this 
dissipation criterion is likely to have more physical mean- 
ing than the criteria used for examp le in classical sand- 
pile an d in more elaborate models fsee lCharbonneau et alJ 
(2001) for a review). However, there is still some doubt 
about the quality of this criterion, as will be discussed 
later (see section l4"T]) . 

Reconfiguration of the field. When the current density ex- 
ceeds a given threshold at some real-space grid points in a 
given plane, the nonlinear terms of equations (11121) which 
are negligible during the loading phase (the off phase) be- 
come large and dominate the dynamics of the fields. They 



are quickly balanced by the dissipative terms when the en- 
ergy cascade reaches the dissipative scale. This "on" phase 
(see section l2~TJ) is modeled by a diffusion- like process for 
the magnetic and velocity fields which tends to reduce the 
magnitude of the current density and the vorticity. 

The detail ed algorithm is an upda ted version of the one 
introduced bv lEinaudi fc Vellil l)l999|) . At a time t for plane 
zq we compute the current density J z (x, y, zq) — — A±A Z , 
where A z = (a + — a~)/2 is the magnetic vector poten- 
tial and Aj_ denotes the Laplacian operator in the plane 
zq. If at some grid point (x,y, Zq) the value of \ J Z \ ex- 
ceeds the threshold J c , A z is updated in the time St c (with 
St c <C St ) according to the equation A z (x, y, zq; t + St c ) = 
A z (x, y, Zq] t) — ri 5t c J z {x,y, ZQ\i), which corresponds to 
current dissipation. The current density J z correspond- 
ing to A z (t + 5t c ) is then computed, and this dissipation 
process is iterated until J z does not exceed the threshold 
anywhere in the plane zq. However, note that after the 
first iteration of this process, we take C J c as a threshold 
instead of J c . The "dissipation efficiency" C is a number 
between and 1 which guarantees that the system is in a 
relaxed state after the whole dissipation process. 

Energy release During this relaxation process, magnetic 
and kinetic energies arc released. The energy release in 
each plane can easily be computed from the variations 
of (£r]_) and (v^_) in the plane. It is the primary variable 
for our statistics. Note that topological modifications of 
the magnetic field may be expected : the connectivity of 
the magnetic field lines is modified because of the field 
diffusion. One of the possible interpretations of this phe- 
nomenon is magnetic reconnection (see however the dis- 
cussion in section |4~TJ> . 

2.3. Time and space scales. 

Let 5x and Sz be the distance between grid points in the 
x (or y) and z directions respectively, and let St be the 
time step. If we assume that the loop length L is 1 to 
100 Mm, then 5z is 30 km to 3 Mm for a typical resolution 
of Nl — 30 points (i.e. 30 planes along z). The analog 
assumptions for a loop width £ (— L/10) of 0.1 to 10 Mm 
give Sx between 1.5 and 150 km for a typical resolution of 
N e = 64. 

We can also determine time scales for the model : as 
the Alfven speed is one in the model units, i.e. in units 
of Sz/St, we have St = Sz/va- the time step is the time 
needed by the Alfven wave to propagate from one plane 
to its neighbors. For B — 10~ 3 to 10 -2 T (i.e. va ~ 1 to 
lOMms -1 with density po « 10~ 12 kgm~ 3 ), this yields St 
between 3 x 10 -3 s and 3 s. Another time scale in the sys- 
tem is the coherence time of photospheric loading Sti. It is 
modeled by a periodic re-initialization of the coefficients 
of the loading increment , which occurs every 200 time 
steps St, or 0.6 to 600 s. This is small compared to obser- 
vational evaluations of the photospheric coherence time, 
but the relevant point is the good separation between time 
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scales ; besides, larger values of the photo-spheric loading- 
coherence time do not alter the results of the model. When 
no other indication is given, the time step St is the unit 
of time ; for example, on Fig. [21 the x-axis range maxi- 
mum is 10 000 St, i.e. between approximately 30 s and 8 h. 
At last, the shortest time scale in the model is the cas- 
cade time scale St c , which is the time step for dissipations 
within a cascade, and which is analogous to the non-linear 
time scale of MHD models. It is assumed to be completely 
separated from the other time scales, i.e. St c <C St <C Sti. 

A direct consequence of the separation between the 
cascade time scale St c and the time St of wave propagation 
between planes is that the helds of neighboring planes are 
expected to be uncorrelated. 

3. Numerical results of the model 

3.1. Model behavior 

The simulations presented in this paper have been per- 
formed on a local quadri-RS/6000 IBM workstation at 
IAS. A typical run of 200 000 time steps with a resolution 
of Nl x Nf = 30 x 64 2 takes between 2 days and 2 weeks 
for one CPU, depending on the parameters. 

Initial growths of energy and dissipation. As the initial fields 
in the simulation box are zero, the initial kinetic and mag- 
netic energies are zero. The loading phase inputs energy 
into the system at each time step St which gives a growth 
of the total energy of the system as shown on Fig. [21 
Then the current density threshold can be reached at some 
points and dissipation occurs, which slows down the ini- 
tial energy growth. At the same time, the average rate of 
dissipation increases until a stationary state is reached. 

Histograms and fitting methodology. The heights of the 
bars of the histograms we plot are normalized by their 
width and they are are divided by the number of events, 
i.e. our histograms are empirical PDFs. A least-squares 
linear fit is then performed in bi-logarithmic axes, on a 
range determined by visual inspection of the histogram 
(see Figure ^p) . This gives error bars on the slope of the 
linear fit, which is the slope of the expected histogram 
power-law. However, we should keep in mind that the 
choice of the fitting range often introduces much larger 
error bars (typically ±0.1 to 0.2) than the error bars of 
the least-squares linear fit of the slope (typically ±0.01 to 
0.05). The error bars we give from now are conservative es- 
timates taking into account the fitting range uncertainty. 

Choice of the variable used for the statis tics and general 
shape of the PDF. Former studies jAletti et all (|200(t : 
lAlettil ll200l|n plotted the histograms or PDFs of the mag- 
netic energy dissipation AE'tot calculated in the whole sim- 
ulation box (Fig. [2J. The global shape of the PDF was a 
Gaussian. A power-law Pr(A-Etot) °c AE^ seemed to ap- 
pear as a deviation from the Gaussianity in the tail of 




2000 4000 6000 8000 10000 2000 4000 6000 8000 1 0000 

t t 



Fig. 2. Initial growth of energy dissipations for parame- 
ters (a) (see Table H): magnetic energies dissipated in the 
whole simulation box (a and b) and in one given plane 
(c and d) are plotted for the 10 000 first time steps of 
the simulation, a) and c) have linear coordinates, b) and 
d) have semi-logarithmic coordinates. Note that no mag- 
netic energy is dissipated in the box until t = 90, and in 
the plane until t = 240. Note also that a stationary state 
is reached from t w 5000 (bottom right). 

the distribution, but it only spanned half a decade, which 
makes it perhaps not so relevant. On the contrary, we 
choose to plot the PDF of the magnetic energy AEi dis- 
sipated in one given plane i £ [1,Nl] (Fig. Ob). As the 
computations are done in the Fourier space, this is our 
primary variable. The power-law that can be fitted to the 
PDF of this variable has a much wider range (more than 
2 decades) and is much less steep (the index is between 1 
and 2) than in the former case. 

This can be explained by a Central Limit Theorem af- 
ter remarking that AI? to t = J2i AEi and that the Ei 's 
are quasi-independent (as expected, the correlation be- 
tween fields in neighboring planes is very low), thus the 
PDF of AE tot is the convolution of the PDFs of all AE t 
for i £ [1,Nl\. The difference between the PDFs in both 
cases stresses the importance of the choice of the variable 
used for the statistics. It also emphasizes that in the case 
of the statistics of observational data, we have to be care- 
ful about the definition of an "event" . 

Both distributions of AEi and of AE show a maxi- 
mum. In the case of AEi, it is a consequence of the fi- 
nite range of the power-law distribution ; the position of 
the maximum depends the average event size and on the 
slope. In the case of AE, knowing the distributions AEi 
in all planes i, it can be seen as a simple consequence of 
the Central Limit Theorem. 
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Fig. 3. Histograms of magnetic energy dissipations a) 
in the whole simulation box ; b) in one given plane. The 
dotted lines correspond to one event per histogram bar 
(0.02 decades wide each). 



Effect of initial growth on statistics. During the initial en- 
ergy growth, the PDF of the energy of events is different 
than during the stationary state. In particular, it is shifted 
to the left, i.e. the events are smaller. As a result, the left 
part of the PDF of events energy gets higher than what 
it would be if stationary state events only were taken into 
account, as seen on Fig. As we are interested in sta- 
tionary state events, we choose to exclude events occuring 
during the initial energy growth from the statistics. 

Typical fields. As the model is built on phcnomcnological 
evolution rules, it is not expected that the fields produced 
by the simulation coincide with any real picture. However, 
as we have tried to be as close as possible to the original 
MHD equations it is interesting to see how far the fields 
are realistic and what the limits of this phenomenologi- 
cal model could be. Typical magnetic and current density 
fields for a = 2 and a = 4 are shown on Fig. [SJ On both 
samples but especially for high values of a, we can notice 
that high current densities occur in magnetic "islands" 
and in regions where magnetic field densities are high. We 
do not observe many structures such as current sheets or 
possible reconnection sites, as will be discussed in section 
14.11 although they are more present for small values of a. 
Large-scale photospheric forcing (large a) leads to large- 
scale structures. 




i<r~ io 5 io 6 10 7 



^^mag tot 

Fig. 4. Effect of initial energy growth on the statistics 
of magnetic energy dissipations in the entire box : one his- 
togram (solid line) only takes into account stationary state 
events, whereas the other one (dashed line) also takes into 
account the events produced during the initial phase. 




10 20 30 40 50 60 10 20 30 40 50 60 




Fig. 5. Top : magnetic and current density fields pro- 
duced in a plane of the simulation box, for parameters 
sets (i) (a = 2, left) and (o) (a = 4, right). On both 
samples, the magnetic field lines are superimposed on a 
grayscale map of J| with large values in black. Bottom: 
magnetic field (left) and current density contours (right) 
issued fr om numerical integratio n of RMHD equations 
( courtesy iGeoreoulis et alJ {n)98|)). 



3.2. Parametric study of event energy PDFs 

A parametric study is performed in order to explore the 
influence of the simulation parameters on the magnetic 
energy dissipations PDFs. A reference set of parameters, 
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Table 1. Sets of parameters used for the parametric study. 
C is the dissipation efficiency, 77 is the magnetic resistivity 
and a is the index of the ID power-law loading spectrum. 
Parameters which are different from parameters set (a) 
are shown in italic font. 
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called (a), is chosen (see Tabled, and it gives the PDF 
shown on Fig. [31 - The PDFs obtained for other sets of pa- 
rameters will be compared to the PDF obtained for (a). 
Most of the sets of parameters correspond to the modifi- 
cation with respect to (a) of one parameter (dissipation 
efficiency C, magnetic resistivity 77, loading spectrum in- 
dex), which is in italic in Tabled All simulations were 
performed on 200 000 timesteps, which seems sufficient 
to achieve a long stationary state after the initial energy 
growth phase, and to achieve good statistics for the PDFs. 
Histograms were done with data from the 100 000 last 
timesteps. 

Other parameters, which are not changed during the 
study, include the grid size (see above) and the current 
density threshold J c = 300. A higher grid resolution would 
have been interesting so as to get a broader wavelength 
range, but it would need a rescaling of other parameters 
and longer computation times. The current density thresh- 
old fixes the scale of current density, so its value has no 
intrinsic meaning. 

Dissipation efficiency. Dissipation efficiency tells how much 
the system gets relaxed after a series of iterative dissipa- 
tions : the current density threshold J c is replaced by a 
new threshold C ■ J c after the first dissipation. With re- 
spect to the value C — 0.5 used in parameters set (a), the 
dissipation efficiency can be set to almost any value of its 
range [0, 1] of valid values with almost no visible change 
in the PDFs. 

Magnetic resistivity. Magnetic resistivity rj could vary in 
the range [3 • 10~ 4 , 3 • 10~ 3 ]. A numerical stability analy- 
sis, given the time step fixed to one and the wavenumber 



range, shows that higher values of rj would result in nu- 
merical instability. On the other hand, lower values of rj 
would lead to longer computational time. However, rj has 
mainly an influence only on the dissipation process length; 
a change in the value of r\ has little influence on the his- 
tograms of dissipated energies. 

Loading spectrum. The reference parameters set (a) has 
a loading spectrum index a = 5/3, corresponding to the 
spectrum in the inertial range of Kolmogorov turbulence. 
In parameters sets (a) and (h) to (o), a varies from 3/2 to 4 
by a maximal interval of 1/3. Another series of simulations 
was performed with lower loading power values to explore 
the influence of the ratio P\ oa ,d/ J c - For both series, for high 
values of a, the power-law is well defined (3 to 4 orders of 
magnitude wide). Its slope index is approximately —1.6, 
and this value depends neither on the loading spectrum 
index a (as seen in Table and Figure El and as will be 
discussed in section l4~T|l nor on the loading intensity -Pi oa d- 
For low values of a, however, power-laws were difficult to 
obtain, and their slopes were sensitive to both loading 
spectrum index and loading intensity. 

Table 2. Variability of the event energy PDF power-law 
index £ as a function of the loading spectrum index a. 
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3.3. Statistics of durations of events 

Durations of events extend over two decades. They are 
indeed a discrete variable, which is a multiple of the cas- 
cade time step 5t c , and their maximum value is a few 
hundreds times St c . Histograms can be obtained, although 
their width is too narrow to perform relevant power-law 
fitting (see figured). 

The duration of an event is correlated with its energy, 
like dEi oc dt}^ 6 , as seen on the scatter-plot on figure |S1 
Another way to visualize this correlation is to select events 
according to their duration, and to plot the histogram of 
energies of events from this population, as shown on figure 
HO One possible observational consequence could be that 
missing long-duration events, due for example to short ob- 
servation times, can lead to energy histograms with nar- 
rower ranges and steeper slopes. 
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Fig. 6. Event energy PDF power-law index £ as a function 
of the loading energy spectrum index a (Table On the 
left, the big error bars are a consequence of the fact that 
the fitting range is not well defined. 



2.16+/- 0.13 




Fig. 7. Histogram of events durations for parameters set 
<1>- 



Fig. 8. Correlation between events duration and energy 
for parameters set (1). 



1.62+/- 0.01 




dE. 



Fig. 9. Event energy histograms for different event dura- 
tions ranges, which form a partition of the whole dura- 
tion range, from low duration (left, dashed) to long du- 
ration (right, dashed). The sum of these histograms (i.e. 
the event energy histogram of all events as in figure 03) is 
shown as a solid line. 



4. Discussion and Conclusions 

4.1. Properties of the model 

The CA model presented here differs from previously pre- 
sented work in two features. On the one hand, the energy 
pumping due to photospheric motion is known quite accu- 
rately via solution of the Alfven wave propagation equa- 
tion with reflecting boundary conditions. Each cell there- 
fore receives and sends energy to neighboring cells along 
the loop axis via a wave equation. Energy redistribution 
to cells on the same loop plane however occurs using an 
instability and redistribution criterion. This criterion is a 
threshold in current, and satisfies the basic requirements 
for magnetic fields of divergence-free conditions and real- 
istic current redistribution. 



slope of these power-laws is quite the same and takes a 
value close to — ( = —1.6. If this behavior is confirmed 
by further parametric studies, it could be interpreted as 
the "unive rsal" behavior of SOC systems, described for 
example in iBak et al. 

It is howeve r inte resting to note that 
iGeorgoulis fc Vlahosl l)l998j) observe in their model 
a variability of the slope of event energy (and also peak 
luminosity and duration) as a function of the loading : 
the slope is steeper when the amplitude distribution of 
the loading increments has a steeper power-law. However, 
this discrepancy is due to the differences in the way the 
system is driven : their loading increments are discrete in 
space whereas our loading has a varying spatial spectrum 
and no power-law amplitude distribution variation. 



Power-law slope of energy distributions as a function of pa- 
rameters. When parameters sets lead to wide and robust 
power-law statistics of events energies, it seems that the 



Reconnection, and the quality of the dissipation criterion. 

Snapshots of the superposition of the magnetic lines with 
the current densities (see figure EJ reveal that the sites 
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of dissipation often do not correspond to reconnection of 
magnetic field lines ; the reason is that one may have in- 
tense currents (and our dissipation criterion is based pre- 
cisely on current intensity) without having a topology of 
the field where reconnection occurs (instabilities such as 
resistive kinks and/or tearing modes are triggered by a 
combination of currents and current gradients). Though 
the typical field structures observed resemble fields from 
turbulence simulations, one sees that our simplified model 
distributes dissipation in a different, more homogeneous 
way. Of course, a cellular automaton model, which mod- 
els the non-linear terms of the equations through sim- 
ple threshold dynamics, is not supposed to generate such 
structures, but we can try to understand what can be done 
to improve the model. The physical quantities available in 
our model make it for example possible to use a more 
elaborate dissipation criterion which would model more 
accurately the reconnection instability threshold, such as 
for example introducing a combined criterion on current 
and current gradients as a trigger for relaxation. 



of observations (~ 100 km at best) compared to the dissi- 
pative scales of the system (« 100m). It is clearly shown in 
this paper that the plot of the PDF of the magnetic energy 
dissipated in a given plane of the model has a shape quite 
different from the PDF of the same observable but for the 
whole simulation box. The bias here consists of a steep- 
ening of the power-law slope with an index greater than 
2. This result suggests that the limited instrumental res- 
olution may be a source o f error as we ll . The SOHO/EIT 
observations analyzed bv lAletti et alJ l|200f)() seem to il- 
lustrate quite well this interpretation : the pixel intensity 
distribution power-laws are steeper (with indices going up 
to 5.6) at lower resolution. Most of the coronal structures 
in the quiet Sun are indeed smaller than the spatial res- 
olution of EIT. Besides, the domain where the power-law 
is fitted on the PDF is reduced which leads to larger error 
on the value of the index. Reliable statistical results could 
be accessible with higher instrumental resolution but also 
by using mathematical tools like for example Pearson's 
method l|Podladchikoval 12002^ . 



4.2. Comparison with observations 

Impulsive coronal events are statistically distributed over 
an energy range of some eight orders of magnitude. Since 
Parker's idea of the existence of nanoflares, it has been 
thought that at some point in the quest of small scale coro- 
nal structures we will break the self-similarity of the solar 
corona by the observation of a steepening of the flare dis- 
tribution with finally a power-law index C greater than the 
critic a l value of 2. Recent d ata analyses (IKrucker fc Benz , 
Il998l lAschwanden et all l200(j IParnel fc JuppI 12000) 
seem to show this behavior with observed values of £ go- 
ing up to 2.6. Empirical formulas have been used here to 
determine flares energies f rom observed luminosities but 
new analyses of the data l|Aschwanden fc Charbonneail 
l2002t iBenz fc Kruckeil h()(m reveal in fact the existence 
of a bias due to the finite range of temperature on 
which the observations are made. The correction of this 
temperature bias leads eventually to a value of £ close 
to 1.6 valid for the whole range of energies, from un- 
resolved, X-ray observations with the Solar Maximum 
Mission (SMM) to Extreme Ultra- Violet observations 
with the Extreme-Ultraviolet Imaging Telescope on the 
Solar and Heliospheric Observatory (SoHO/EIT) and the 
Transition Region And C oronal Explorer (TRACE) - see 
lAschwanden et all l|200ffl . Such an observational bias, as 
well as the one described in the end of paragraph ia. Al show 
the importance of defining well what an event is and what 
its characteristics are, when observing the corona but also 
when using statistical flare models. 

The comparison between model distributions and ob- 
served distributions is actually not an easy task, even if it 
is tempting to compare the power-law slope 1.6 obtained 
by our model to the observed global slope of 1.6. The first 
pitfall for comparison between statistical results of obser- 
vations and models may be linked to the spatial resolution 



However, we think that resolution has not always such 
a dramatic effect on the slope of the PDFs, and that it 
is still interesting to model statistics of individual events 
(i.e. in one plane, in the case of our model). The con- 
vergence to a Gaussian when summing the energies of 
events before doing statistics, predicted by an argument 
lying on the Central Limit Theorem, may indeed be much 
slower in the case of observed, real micro-events than in 
the ideal case of independent events with low-moments 
distributions: the distributions of real events energies are 
much wider than the modeled distributions and their mo- 
ments may be greater. As a result, depending on the ob- 
servational conditions, there could still exist a quite wide 
power-law after summation (due to lack of resolution) of 
a large but not too large number of events, and the slope 
of this distribution may be still close to the slope of the 
original distribution of individual events. In this context 
the slope £ = 1.6 we obtain is rather encouraging. 

Another prediction of our model is that durations and 
energies of events scale like dEi cx dt} 76 , or, equivalently, 
dti cx dEl f1 - 76 = d E?- 57 . This exponent 0.5 7 is in quite 
good agreement with lBerghmans et alJ l)l998|) . who report 
that observed events durations scale like their radiative 
loss at the power 0.5. 

We should however emphasize that statistical flare 
models usually give energy dissipations, while observa- 
tions give luminosities at some given wavelengths and the 
infered energies depend on models. It is therefore crucial 
in the future to develop models including the produc- 
tion of observable quantities, in order to provide stronger 
links between models and observations but also to quan- 
tify more precisely the weight of observational biases. The 
first agreements obtained during the last decade between 
statistical predictions made by theoretical models and ob- 
servations are however very promising. 
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4.3. Summary and conclusion 

In this paper we have presented a three-dimensional sim- 
plified model inspired by th e RMHD equatio n s who se first 
version was introduced bv lEinaudi fe Vellil l)l999lh This 
model mimics a coronal magnetic loop anchored in the 
photosphere whose footpoints are driven randomly by con- 
vective motions. The slow driving of the magnetic foot- 
points leads to storage of energy along the coronal loop 
and eventually to dissipation through impulsive events. 
The characteristics of the model are the following : (i) the 
model describes Alfven wave propagation along a loop ex- 
actly ; the internal structure of the loop is described by a 
set of planes distributed along the axis and ending in the 
photosphere from which the information propagates ; (ii) 
the external forcing applied to the two boundary planes is 
expressed as a turbulent spectrum in Fourier space ; (iii) 
when the criterion of instability is satisfied, the dissipation 
of current density and vorticity takes place non-locally in 
Fourier space but still locally in physical space ; (iv) Fast 
Fourier Transforms (FFT) are implemented in the numer- 
ical code to use the dual physical/Fourier space. 

A numerical study has allowed to quantify the role 
of the parameters, especially forcing, in the behavior of 
this model and in statistical properties of coronal events. 
The slope of event energy histograms was found to be al- 
most constant, in accordance with a SOC-like "universal" 
behavior, and consistent with the values given by obser- 
vations. Event durations statistics were performed, and 
correlations with event energies are also compatible with 
observations. Different possible observational biases were 
pointed out, all of them resulting in a narrower power- 
law range on histograms and in a steeper slope than in 
the statistics of all elementary events : a bias due to the 
limited spatial resolution, which gives a possible interpre- 
tation of recent observations made with the instrument 
EIT on board SoHO, and a bias due to limited observa- 
tion durations. 

Improvements are always possible to provide a better 
description of the physics. One can imagine some ad-hoc 
rules to obtain for example a correct picture of the recon- 
nection process. But the probably most interesting (and 
difficult) study is about the incorporation of the non- linear 
dynamics in a more realistic way than the simple on-off 
mechanism of CA models, but without going directly to 
the MHD equations. From a pure observational point of 
view it seems crucial to have as precise as possible an esti- 
mate of the possible biases to determine the effective value 
of the power-law index £ of the energy distribution. Indeed 
the confirmation of the sub-critical value of £ could be a 
serious challenge to Parker's hypothesis of coronal heating 
by a swarm of nanoflares. 
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